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分 数 阶 Nernst-Planck 方 程 的 有 限 差 分 / 谱 元 法 求解 * 


EWB, 许 传 炬 
(厦门 大 学 数学 科学 学 院 ， 厦 门 361005) 

摘 ”要 : Nernst-Planck 方程 是 用 来 描述 在 离子 浓度 梯度 VC 及 电场 VV 共同 存在 的 情况 下 ， 穿 过 渗透 膜 
的 离子 (如 钙 ， 钊 ， 钠 ， 氯 ， 镁 等 ) 流 了 的 方程 。 但 是 ， 计 算 Nernst-Planck 方程 的 数值 解 会 遇 
到 一 些 困难 。 本 文 考虑 用 以 描述 神经 细胞 中 离子 反常 扩散 现象 的 电缆 型 简化 的 分 数 阶 Nernst- 
Planck 方程， 提出 了 一 个 时 间 有 限 差分 /空间 谱 元 法 对 该 方程 进行 数值 求解 。 我 们 给 出 了 数值 方 
法 的 详细 构造 过 程 以 及 实现 方法 。 结 果 表 明 数 值 解 在 空间 方向 上 具有 指数 阶 收 敛 精度 ， 在 时 间 
方向 上 具有 2 一 a 阶 精度 。 最 后 ， 通 过 计算 一 个 具有 实际 背景 参数 的 问题 说 明 所 提 方 法 的 潜在 应 
用 。 
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细胞 生命 活动 过 程 中 伴随 的 电 现 象 ， 主 要 由 存在 于 细胞 膜 两 侧 的 电位 差 ( 称 为 膜 电 
位 ，membrane potential) 来 刻画 。 跨 膜 电 位 差 形成 的 一 个 主要 原因 就 是 细胞 膜 上 存在 选择 
性 离子 通道 ， 大 多 只 允许 特定 的 离子 通过 。 在 电 扩 散 过 程 中 ， 一 方面 ， 根 据 Fick 定律 离子 要 从 
浓度 高 的 区 域 往 浓度 低 的 区 域 扩散 ; 另 一 方面 ， 由 Kohlrausch 定律 ， 带 电 粒 子 在 电场 力 的 作 
用 下 有 漂移 速度 。 在 这 两 种 力 的 作用 下 ， 神 经 细胞 中 离子 的 流通 量 为 由 


Ji = -u (VC + aF avv), 


其 中 入 是 电势 养 ， 拓 为 离子 i 的 流量 ， 本 表示 扩散 常数 ，Ci 是 离子 浓度 ， 竺 是 化 合 价 ， 环 是 法 
拉 第 常数 ， 忆 是 理想 气体 常数 ， 工 是 绝对 温度 。 
每 一 种 离子 都 遵守 质量 守恒 定律 ， 于 是 得 到 下 面 的 Nernst-Planck 方 程 (N-P 方程 ) 
OG 


Be + Y=0, 


实际 的 Nernst-Planck 方 程 是 三 维 的 ， 计 算 难 度 较 大 。 为 了 简化 模型 ， 许 多 研究 考虑 细 长 神经 
元 情形 ， 引 进 所 谓 的 电缆 (Cable) 模型 。 如 Qian 和 Sejnowskil 所 提 到 的 ， 为 了 使 问题 简化 ， 
我 们 考虑 一 个 细 长 的 电缆 模型 ， 假 定 轴 方 向 的 电流 和 离子 浓度 在 圆柱 的 横 截 面 上 是 一 致 的 ， 并 
且 横 截面 方向 的 电流 只 在 圆柱 表面 存在 与 绕 圆柱 轴 心 的 辐 角 无 关 。 这 些 假设 使 得 电 扩 散 问题 简 
化 为 没 着 圆柱 轴 心 的 一 维 问题 ， 此 时 在 刺激 电流 J,; 下 方程 可 写 为 站 
Oc; OC, zuviF 0 OV 4 
Gat Rr Ba Cae ) ra 
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EP Imi 是 离子 i 的 电流 密度 (规定 向 外 为 正 )，1s;i 是 刺激 电流 密度 ，d 是 电线 圆 截面 的 直径 。 
为 了 使 系统 封闭 ， 要 加 上 一 个 与 电势 差 相 关 的 方程 
V(r, 9=W+ 2 Ce (i (2) 
EF vV 是 初始 电势 差 ，CGCio 是 离子 i 的 初始 浓度 ，Cm 是 单位 面积 的 细胞 膜 电 容 。 
方程 (1)-(2) 的 计算 已 有 许多 成 熟 的 方法 。 最 近 ，Barkai 等 人 四 根据 实验 中 发 现 的 反常 扩散 
现象 导出 分 数 阶 的 Fokker-Planck 方 程 。 事实 上 ， 更 早 一 些 的 研究 如 Metzler 等 人 里 也 发 现存 
在 类 似 的 现象 。 初 边 值 问 题 的 分 数 阶 Nernst-Planck 方 程 具有 如 下 形式 


ac; _ Rpil-a 8? Ci znF ð ov _ =D Pa "i 
Ot = D; (ug De 2 RT z(a =) Faz) Isis (3) 
且 满足 初始 条 件 
Ci(z,0) =Cio(a), V(x,0)=Vola), VreQ, (4) 
和 边界 条 件 
Ci(z,blano = Cio(z)lan, Y(z,blen = Vo(z)len, VEE T, (5) 


其 中 RDI, 0 < a < 1， 表 示 关 于 时 间 变 量 上 的 1 — aft Riemann-Liouville Rot FR) 


1 d ff vdr 


rh En vt € [0,T). 


aD} -ev(t) = 


我 们 考虑 K+ 和 Nat 两 种 离子 的 情形 ， 即 i = K+，Na+。 下 面 描述 的 方法 可 以 直接 应 用 到 
多 种 离子 情形 。 这 两 种 离子 的 膜 通 量 如 下 [9] 


Im, Nat = Inam h(V = Via), Im K+ = gen (V = Vk), 


SEP gy, 和 gi 分 别 表示 最 大 钊 和 钠 电 导 , Vna AV, 分 别 是 钠 和 钾 平 衡 电势 。 无 量 纲 变量 m, n, h 
满足 如 下 的 一 阶 常 微分 方程 


dm 


dn dh 
“dt. ry = Qn(1 =n) — Grn, dt 一 ah 人 (1 ~h) — Pah, 


= Qm(1 —m) ~ bmm, q 


其 中 


0.1(V + 25) _ 0.01(V + 10) 


a a 0.07exp (元 )， 
SS + Pers Ty, ee | h = 
exp( 4125 — 1) exp(+#/2 — 1) 20 


Qm = 


V V 1 
Bm = 4exp (=): Bn = 0. 125 exp (+=), Gn = exp( 240 +1) £30 44) 


表 1 列 出 了 各 种 量 纲 和 求解 方程 所 需 的 一 些 参数 。 初 值 条 件 是 静止 状态 下 的 条 件 。 
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Rl: N-P 方 程 中 未 知 量 的 量 纲 和 一 些 参 数 


G 9.6485 x 10’mc/mol 
Vya 3.47 x 107cm? /ms 0.005ya 
x 5.11 x 10-8cm?/ms 0.005 pa 
2 A mol 
A 


将 Riemann-Liouville 44 MAA AT 1 A 作用 于 (3) 的 两 边 ， 并 利用 分 数 阶 微 积分 组 合 性 
JR By!) 


2C zuF ð DY 4 as 
PC — — Seal Gas) = apr mi li i=Kt, Nat, 
Or > RT 5a (Ce ) dFa ™i ‘sh a ó 

6 


Fd 
V=Wi+t 1G XG = Ci,0) Zi; 


其 中 De 表示 a 阶 Caputo 导数 加， 无 ; = Ie Isio 

目前 已 有 一 些 关 于 整数 阶 Nernst-Planck 方 程 数值 计算 的 研究 工作 。 如 Qian? 提出 时 间 方 
向 用 一 阶 差分 ， 空 间 方 向 用 中 心 差分 对 Nernst-Planck 方 程 进行 离散 。Samson 等 人 [中 用 有 限 差 
分 /有 限 元 法 求解 扩展 的 Poisson-Nernst-Planck 方 程 。 但 是 分 数 阶 方程 的 求解 与 整数 阶 相 比 有 具 
有 本 质 的 不 同 。 本 文 将 提出 一 个 在 时 间 方 向 上 采用 高 阶 有 限 差分 法 ， 空 间 方 向 上 采用 谱 元 法 对 
分 数 阶 Nernst-Planck 方程 数值 求解 。 


2 ”有限 差 分 / 谱 元 法 


本 节 详 细 描 述 我 们 的 方法 。 
2.1 时 间 离 散 : 有 限 差分 格式 
对 于 给 定 的 正 整数 M > 0， 令 如 = nAt, n= 0,1,… ,M， 其 中 时 间 步 长 At = 2. 那么 对 
所 有 的 0 < n < M， 我 们 采用 如 下 差分 格式 逼近 分 数 阶 时 间 导 数 
1 z ee 0, f(T) 
| 


A = Fa) Ooh, Gen 7 


1 
r-a) F 


£ ; tj+1 
人 Hi) J (tn41 — 7) -dr + Qt) 
tj 


n 
=0 


210 I E 数 学 学 报 第 27 卷 


= Rep Ta Le (tees) — HG) (+ 1-3 aae) ra 
j=0 


= KT a YO (ji 一 Fn)(G1L+ 用 = 94) H rat 
j=0 


Bo >. bs ( F(te—341) — fltn-;)) +7873, (7) 


j=0 


i 


这 里 bj) = (j +1) — jte, Bo = ATOA rn+1l 为 截断 误差 项 。 
引 理 1B@ 格式 (7) 的 截断 误差 为 


rati < OAt >). 


将 差分 格式 (7) 用 于 方程 (6) 中 的 DEC(a, tri) 项 ， 将 二 阶 Adams-Bashforth 格 式 用 于 Ini 
项 ， 其 它 的 非 线 性 项 进行 隐 式 处 理 ， 得 到 如 下 的 半 离 散 格式 


n+1 一 1 i 
port — ZA _ BHF Ê (opm VO) A Po 0s =t 


Ox? RT Oz Ox 
— Bobn V? + ee OE -pP -ẸP, i=Kt,Nat, (8) 
Fd 
yor oa you PTa > g CO0)zi， 


Soh V4} BIE V (a, tn41), CPE MIE Ci(x, tn41), V° = Vo, CP = Cho。 方 程 (8) 以 及 边界 条 
件 (5) 和 初始 条 件 (4) 构成 一 个 完整 的 半 离 散 问题 。 

2.2 ”空间 谱 离散 

为 了 更 好 地 模拟 在 给 定 的 位 置 注入 刺激 电流 引起 的 离子 扩散 ， 我 们 用 Legendre 谱 元 法 进行 
空间 离散 。 使 用 谱 元 法 的 优点 之 一 是 : 通过 区 域 剖 分 ， 可 以 选择 适当 的 剖 分 点 使 得 在 刺激 点 附 
近 有 足够 多 的 节点 来 模拟 刺激 电流 。 

W L(A) 为 A 上 的 Lebesgue 平 方 可 测 函 数 空间 ， 相 应 的 内 积 为 (u,v) = fa uvde. H™(A), 
Hy" (A) 为 标准 的 Sobolev 空间 。 

考虑 问题 (8),(5) 的 弱 形式 : ROCH € (H(A), Veti e 所 (A)， 使 得 对 任意 的 W e€ 
H(A WE 


ocnt | 
Oz Ox 


+ ay (cpt ores, ow) = (yn W), i=K+,Nat, (9) 


Ox 


Bo (Cr+, W) +r ( 


Fd 
vt 了) 一 一” Crtla, W ) =(V°,W), 
( ) TAOD = ) mM) 
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其 中 
n-1 , 4 p 
YH = By X (b; — bj )CP? + BobnV® 一 IFz elimi ~ TRG") ~ Ls 
j=0 


Fd 
0 0 Co A 
ja o> ee 


我 们 将 要 构造 的 谱 元 法 基于 标准 的 Galerkin 方 法 。 为 此 ， 先 进行 网 格 剖 分 ， 在 区 问 A 中 选 
取 太 十 1 个 点 ， 从 左 到 右 依 次 编号 为 a = ao <a, <… < ak = b， 这 样 区 间 人 就 被 分 成 个 
互 不 相交 的 子 区 间 Ap = (ak_1;,ak)， 它 们 满足 天 = UK Ak。 定 义 A = (一 1,1) 为 参考 区 间 ， 显 
然 存 在 一 个 从 A 到 Ax 的 一 一 映射 Fy 


z=F(r)=l(r+1)/2+ag1, reA, lk = ak — ak1. 


Galerkin 谱 离 散 基本 思想 在 于 用 高 阶 多 项 式 来 逼近 解 。 为 此 ， 定 义 A 上 的 分 请 多 项 式 空间 

为 
Py (A) > {v | vo F(A) € Py(A), k= Tee ,K}, 

其 中 Py (A) A Â EAR N 的 多 项 式 空间 。 

ABU A) Vy = P(A) U H(A). BI E H g (9) 的 Galerkin 谱 离散 如 下 ， ROTH E 
V, VET e Vy， 使 得 对 于 任意 的 Wn € Vw 满足 

och BT 
ðr ”pr ) 


Bo (COR, Ww) + n( 


ZF (mii OVE OW. < , 
RT (crx a 3 Fan) = (MN, Wn), i=Kt, Na", (10) 


(EH, Wa) - gee (X Ak a Wu) = (09, Wi), 


其 中 
n < n—j 4 n n— 7 
YN? = Bo 2 6; J bia) CEN” + Bobn V? = IFz Cmi <6 Ina) — Isi. 
J= 
3 ”计算 实施 


Galerkin 383 ja] H (10) 实现 起 来 比较 困难 ， 因 为 精确 计算 (10) 中 的 积分 很 费时 。 一 种 代替 
方案 是 用 数值 积分 取代 精确 积分 。 由 于 Gauss-Lobatto 数值 积分 所 具有 的 高 精度 ， 我 们 选择 用 
它 近 似 (10) 中 的 积分 。 在 区 间 A 上 的 Gauss-Lobatto 积 分 公式 为 


N 
feoir ~ Seu, 


i=0 
其 中 是 Gauss-Lobatto-Legendre (GLL) 点 ，wi; 是 相应 的 积分 权 值 。 上 述 积分 对 任意 的 $ E 
Pon—1(A) 是 精确 成 立 的 。Ak 上 的 GLL 积分 可 通过 变量 变换 得 到 


lk ` lk 
(u,v)A, = [ow o Fr(r) 3 dr X (uv) o Felgi) 3 oi 


7=0 
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wet = Feli), we = wi WE A Ar KI GLL 点 。 定 义 人 A 上 上 的 离散 内 积 为 


K 
(u,v)n = So, V)Ak,N， 


k=1 


其 中 


N 
(u,v) ann = 》 ulek yolk wf. 
i=0 
将 (10) 中 的 内 积 用 离散 的 内 积 代替 ， 得 到 一 个 近似 问题 : ROTH EVR, V E Vv， 使 得 对 
于 任意 的 Wn E Vn» 满足 


OCr'# OW, 
n+l ; i,N N 
Bo( Ci Wr) +™( Or ` Ox E 


zin EF n avnt Ow; n 3 
RT (Crs dn i ar a = (Yin Ww) yw, i= K*,Na’, (11) 
n Fd n [7 
(VEH, WN) 一 IC (O Cea a Wn) y = (V°,Wy) vy : 


由 (11) 可 推导 最 终 求解 的 代数 方程 组 ， 具 体 如 下 定义 


N r= 
hi(r)= JI AAA i=0,1,---,N. 
j=0,j#1 J 


显然 及 满足 hi(&j) = 6ij (65 为 Kronecker 符 号 )， 且 对 任意 的 vw € Pw(A), A 


N 
un(r) = D un (&:)hi(r). 


i=0 


S Fol (a) WF. ARUN, 定义 内 (Zz) = hio Fy (2) ABS, hk 有 如 下 形式 


A — gk 
h(a) = Il : 2A zEAk, i=0,1,---,N, k=1,---,K. 
j=0,jzi SË T$} 


FIE hk (EF) = bi50 RE (x) BRA Vn E Ak 上 的 一 组 局 部 基 。 将 ORE) 和 VR+! 用 这 些 基 函 数 展开 


K N 
Cra) =>" > A Erh (e), 


k=1 i=0 


K N 
Veta) =X >》 VR (Et) ht (a), 


k=1 i=0 


K N 
Ya) =Y JYP (Eth (a), i= K+, Nat, 
k=1 i=0 
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代入 (11)， 并 让 检验 函数 Ww 取 遍 所 有 的 由 (z)) 7 = 1,… ,NN 一 1, 上 == 1,… ,K， 我 们 得 到 一 
个 非 线性 方程 组 


Meci+DTBDcei+DTIGIPDo =y, i=Kt,Nat, 
(12) 
Mv- Qici = v, 


其 中 
Gt (Ot! (0) On (ee On (eh) 3 ON (EN) , TKt, Nat, 
v = (VREHEN VET), VETER) , VR ER), 
是 (N 十 1)K 阶 未 知 量 。w9，y; 类 似 定 义 
v = (VE od VCE), VERN), 
yi = (YERE ,VN Ew ,V1 ER lt)". 
M, Qi, Ei, Gi 为 ((NN 十 1)K) x (N 十 1)K) 阶 对 角 和 矩阵 ， 即 
M = diag(M*), M*= diag(Bowt), k=1,---,K, i=0,---,N, 


; | / Edzi : ; E 
Qi = diag(Qt), QF =diag( Te), A =diag(A'), Př = diag( CPR (Et), 
m 


E; =diag(E*), E¥ =diag(yw*), Gi = diag(G*), Gf = diag( 
万 为 ((N +1)K) x ((N +1)K) 阶 块 对 角 矩 阵 ， 即 
D = diag(D*), Di, = Ochf(€F), k=1,---,K, ij=0, N. 


方程 (12) 是 一 个 非 线性 非 对 称 系统 ， 我 们 用 Jacobian-free Newton-Krylov 方 法 回 求 解 。 为 
了 便于 应 用 该 方法 ， 重 新 写 方程 (12) 如 下 


Fi(u) := Me, + DT E,De, + D?G;P,Dv -y =0, i=Kt,Nat, 


F(u) :二 Mv 一 >》 Qici = v? = 0, 


T 
u= ( cns， CK, v) . 


用 Jacobian-free Newton-Krylov 方 法 解 非 线性 系统 (12) 的 过 程 如 下 。 
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首先 ， 给 定 一 个 初始 值 w( 中 ， 用 Bi-Conjugate-Gradient-stab (Bicgstab)l19 算 法 解 关 于 增 | 
E óu 的 线性 系统 


Km BERGER, Te (wu) 是 函数 五 在 w(t" 的 Jacobian 和 矩阵 。 实 际 计算 中 我 们 不 需要 
形成 矩阵 .Tri (ul). BLE, Bicgstab HA RENWAL ARR. dA (ul, bu) 为 函 
BF UM 点 关于 654 方向 的 Fréchet 导数 ， 则 

Ir (u™) du = dF (u™, du) 
Fi (ul + edu) — Fi(ul™) 


< 一 0 € 


_ M(c™ +eôci) -Me™ DTED(c™ + ebc) — DTE,De!™ 


Ei DTGi(B + 6P,)D(v™ + edv) - DTGiRDu(m) 
OD na 
< 一 0 € 


= Méc; + DTE Déc; + D’ GP, Dév + DTIGiP, Dée;, 
Tr, (u™) du = dF, (ul, du) 


Fy (u™ + edu) — F, (uo) 
< 一 0 E 


(m) (m) 

M(v™ + eôv) ~Myu™ = Qi(e;” + €6C;) a Qic;” 
= lim 一 -一 ln 一 一 
< 一 0 € < 一 0 € 


= Mi 一 >》 Qide;, 


其 中 P, = diag(P*), PF = diag(Vpt} (é*)), k=1,--- ,K, i=0, ,N. 
最 后 ， 我 们 根据 


ul) 一 u™ 4+ wou, 


BEEBE ult), HAFO < w < 1 是 用 来 全 局 化 上 述 迭 代 法 的 收敛 性 ， 最 简单 选择 ww 的 准则 是 
使 非 线性 残 量 减 小 ， 即 
| Fut] < (1 - Aw) || F(u™)|. 


取 入 = 10-!， 我 们 测试 w = 1.0, 1/4,1/16,---, AFW EEREN. Newton-KryloviFit 


终止 条 件 是 
õu 
lu] 


在 计算 中 ， 取 tol = 107, RAHE HAHAE REA RAH DRR, AREE 
为 Bicgstab iK (WEN 4B 


< tol. 
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4 ”数值 实验 
算 例 1 为 了 验证 数值 方法 的 有 效 性 ， 我 们 首先 考虑 已 知 精确 解 的 问题 


OPC; 0 OV) _ . 
DeG FF - =(Az) =f, O<x<3, t>0, i=1,2, 
2 
V= >》 G, 
i=l 
SEP ERR SERERA 


V = tt? sin(nz/3), Ci =t?sin(mz), C2 = t! sin(27z/3). 


首先 ， 我们 考察 空间 方向 上 的 误差 。 为 此 ， 先 固定 At = 0.0001， 使 得 时 间 方 向 上 的 误差 足 
够 小 。 取 a = 0.9, 五 = 1， 我 们 考察 C 和 VV 在 空间 方向 上 的 LC? M H! 范 数 误差 的 衰减 情况 。 
从 图 1(a)-(b) 可 看 出 ， 在 半 log RE F (WAER log), L? 和 五 1 范 数 误差 与 多 项 式 阶 数 N E 
线性 关系 ， 这 表明 数值 解 在 空间 方向 是 按 指数 阶 收敛 于 精确 解 的 。 


t=1.0, a=0.9 t=1.0, a=0.9 
= pe f 
L nom 一 一 


~r 一 一 二 一 一 
1 ' i _ Lf nom 一 一 
H seminorm -+~ H seminorm ~+- 


Al: 误差 随 NN 的 变化 情况 


其 次 ， 我 们 考察 时 间 方 向 上 的 误差 衰减 速度 。 为 此 ， 先 固定 多 项 式 阶 数 N = 30， 单 元 个 
BK = 1， 这 样 空间 方向 上 的 误差 可 以 忽略 不 计 。 分 别 取 wa = 0.6，0.9， 在 图 2(a)-(b) 中 ， 
可 看 到 在 log 尺度 下 (对 横 坐 标 和 纵 坐 标 同时 取 1log) L? 和 五 ! 范 数 误差 与 时 间 步 长 At 时 线性 关 
系 。 这 表明 精确 解 与 数值 解 之 间 的 误差 在 时 间 方 向 上 是 代数 衰减 的 ， 且 是 At2-* 阶 收敛 的 ， 这 
与 预计 结果 相符 合 。 

算 例 2 我 们 计算 具有 实际 背景 参数 的 N-P 方 程 ， 计 算 中 所 用 到 的 参数 见 表 1。 在 实际 计算 
中 ， 我 们 计算 Cio — Ci, Vo 一 V 的 值 ， 从 而 考察 离子 浓度 和 膜 电位 基于 静 息 状态 C = Cio, V = 
Vo 时 的 变化 情况 。 为 了 使 Ci,o 一 Ci, Vo 一 V 满足 齐 次 边界 条 件 ， 在 计算 中 我 们 选取 足够 大 的 计 
算 区 域 Q = (一 4,10)， 并 将 其 平均 剖 分 为 14(K = 14) 个 单元 ， 每 个 单元 用 25(N = 25) 阶 多 项 
DOE, ERME r = 0 点 注入 电流 时 ， 其 位 置 恰好 在 第 四 个 单元 和 第 五 个 单元 交点 处 ， 有 
足够 多 的 点 来 模拟 刺激 电流 天。 
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t=1.0, a=0.6 


t? norm 一 一 


H! seminom 
Line of slope 1.4 ——— 


图 2: 误差 随 人 At 的 变化 情况 


当 a = 1.0 时 ， 分数 阶 N-P 方 程 (6) 简化 为 传统 的 N-P 方 程 (1)-(2)。 我 们 采用 二 阶 的 向 后 差 
分 / 谱 元 法 格式 对 N-P 方程 数值 求解 。 图 3 是 Qian 四 在 时 间 上 用 一 阶 差 分 ， 空 间 上 用 中 心 差分 
对 N-P 方程 数值 求解 的 结果 。 取 时 间 步 长 At = 10-?， 图 4 (a), 4 (e) 描述 了 神经 细胞 中 的 膜 电 
位 在 空间 方向 上 和 时 间 方 向 上 的 传导 。 两 者 吻合 很 好 。 

当 0 < a < 1 时 ， 我 们 采用 本 文 所 提出 的 2 — a 阶 有 限 差分 / 谱 元 法 格式 求解 分 数 阶 N- 
P 方 程 ， 研 究 神 经 细胞 中 离子 反常 扩散 时 膜 电 位 在 时 间 和 空间 方向 上 的 传导 情况 。 由 于 
膜 电位 的 传导 速度 随 a 的 减 小 而 加 快 ， 这 就 意味 着 V"+1 与 V* 的 差 值 随 着 a 的 减 小 而 变 
大 ， 而 我 们 把 好 步 的 计算 结果 WV” 作为 用 Newton 连 代 法 求解 V7?+1 的 初始 值 ， 所 以 为 了 保 
证 Newton 和 迭代 法 的 收敛 速度 ， 当 a 减 小 时 ， 取 较 小 的 时 间 步 长 At = 10-?。 我 们 分 别 取 a = 
0.9, 0.8, 0.7， 图 4(b)-(d) 中 描绘 了 在 x = 0,1,2,3 不 同位 置 膜 电位 随时 间 变 化 的 传导 情况 ， 
以 及 图 4(f)-(h) 中 描绘 了 在 上 = 0.5, 1.0, 2.0, 3.0 不 同时 刻 膜 电位 随 空 间 变 化 的 传导 情况 。 从 
4 (a)-(h) 中 ， 我 们 看 到 神经 细胞 中 离子 反常 扩散 时 膜 电位 的 传导 速度 随 a 的 减 小 而 加 快 。 


Membrone potential (mv) 


图 3: 文献 [2 的 结果 : 不 同位 置 的 膜 电位 随时 间 变 化 情况 
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B4: 膜 电位 的 传导 速度 随 a 的 变化 
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A Finite Difference/Spectral Element Method for the Fractional 
Nernst-Planck Equation 


LI Xian-juan, XU Chuan-ju 


(School of Mathematical Sciences, Xiamen University, Xiamen 361005) 


Abstract: The Nernst-Planck equation describes the flux of ions (for example, calcium, potassium, 
sodium, chloride, and magnesium etc.) through a diffusive membrane under the influence of both 
the ionic concentration gradient VC and electric potential VV. However, numerical approximations 
to the Nernst-Planck equation suffer from several difficulties. In this paper, we first briefly recall 
the derivation of the fractional Nernst-Planck equation in a cable-like geometry, which describes the 
anomalous diffusion in the movement of the ions in a neuronal system. Then a method combining 
finite differences in time and spectral element methods in space is proposed to numerically solve the 
underlying problem. The detailed construction and implementation of the method are presented. Our 
numerical experiences show that the convergence of the proposed method is exponential in space and 
(2 — a)-order in time. Finally, a practical problem with realistic physical parameters is simulated to 
demonstrate the potential applicability of the method. 

Keywords: fractional Nernst-Planck equation; spectral method; finite difference method; Newton- 


Krylov iteration 
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